# Set parameters for abundance distrilim_unif <-100par1 <-c(0,log(10), log(3), 0) # the median is exp(par1)par2 <-c(lim_unif,log(1.5), log(1.5), log(10)) # CV is more or less sqrt(exp(sigma^2)-1)names(par1) <- abund_distrinames(par2) <- abund_distri
Below are the abundance distribution functions:
Code
x <-seq(0, lim_unif, by =0.01)glist <-list()for (t innames(par1)) {if (t =="1_unif") { y <-dunif(x, par1[t], par2[t]) tit <-paste0(t, " (min = ", par1[t], ", max = ", par2[t], ")") } else { tit <-paste0(t, " (meanlog = ", round(par1[t], 2), ", sdlog = ", round(par2[t], 2), ")") y <-dlnorm(x, meanlog = par1[t], sdlog = par2[t]) } df <-data.frame(x = x, y = y) g <-ggplot(df) +geom_line(aes(x = x, y = y)) +ggtitle(tit) +theme_linedraw() glist[[t]] <- g}wrap_plots(glist)
# Select the experiments for which delta should varyexp_vary_delta <-"ninter"
Code
# Create a dataframe of legends for the plotslegend_df <-data.frame(exp =unique(paramdf$exp))legend_df$legend <-c("Number of interactions","Species abundance distribution","Mean consumers niche breadth","Standard deviation of consumers niche breadth")
Simulations
perform_simu is FALSE so the simulations will be loaded from the files simres_delta0.rds, simres_delta0.2.rds, simres_delta1.rds.
Code
simres_list <-list()for (i in1:length(delta_list)) { delta <- delta_list[i]print(paste("Delta", delta, "----------------")) delta_name <- delta_names[i]if (perform_simu) { simres <-list("cor"=list("mean_cor"=data.frame(),"sd_cor"=data.frame(),"sd_area_cor"=data.frame()),"example"=list()) rtime <-system.time({if (delta == delta_def) {# Do all simulations jchoice <-1:nrow(paramdf) } else {# Do only the ones in delta_vary jchoice <-which(paramdf$exp %in% exp_vary_delta) }for(j in jchoice) {# Simulation parameters --- parms <- paramdf[j, ]print(paste0("Experiment ", parms$exp, " (", parms$id, ") ---"))for(repi in1:nrep) {print(paste("Rep", repi))# parms$rep <- repi # For debugging# Simulate data ---# Generate consumer and resource abundances # These abundances must be unsorted, because in the simulation# algorithm we sort the resource following their traits and if the abundances# are sorted as well it creates a correlationif (parms$abund_distri =="1_unif") { consumerab <-runif(parms$nconsumer, par1[parms$abund_distri], par2[parms$abund_distri]) resourceab <-runif(parms$nresource, par1[parms$abund_distri], par2[parms$abund_distri]) } else { consumerab <-rlnorm(parms$nconsumer, meanlog = par1[parms$abund_distri], sdlog = par2[parms$abund_distri]) resourceab <-rlnorm(parms$nresource, meanlog = par1[parms$abund_distri], sdlog = par2[parms$abund_distri]) }# Simulate data sim <-compas2d(nconsumer = parms$nconsumer, nresource = parms$nresource,le_grad = le_grad, col_prefix ="b",row_prefix ="p",rowvar_prefix ="tp",remove_zeroes =TRUE,consumer_ab = consumerab,resource_ab = resourceab,ninter = parms$ninter,ratio_grad = ratio_grad,mean_tol = parms$mean_tol, sd_tol = parms$sd_tol, delta = delta,buffer = buffer)# Format simulated data --- comm <-data.frame(sim$comm) consumer_niche <- sim$colvar resource_traits <-data.frame(sim$rowvar)# Multivariate analyses ---# CA ca <-dudi.coa(comm, scannf =FALSE, nf =min(dim(comm)))# Reciprocal scaling rec <-reciprocal.coa(ca)# Measure correlations between mean positions --- res_niches <-get_niches(rec, consumer_niche = consumer_niche, resource_traits = resource_traits,comm = comm,rowname ="resources", colname ="consumers") res <-get_cor(res_niches) ressim <-add_params(res, parms)# Get eigenvalues neig <-10# We fix the number of eigevalues to 10, in case there are less eigenvalues there are NAs in the columns ca_eig <-as.data.frame(t(ca$eig[1:neig]))colnames(ca_eig) <-paste0("l", 1:neig) ca_eig[, names(parms)] <- parms# Chi-squared test sum_eig <-sum(ca$eig, na.rm =TRUE) n_comm <-sum(comm) testval <-chisq.test(comm, simulate.p.value =TRUE) pval <- testval$p.value# Add these to the results ca_eig$sum_eig <- sum_eig ca_eig$n_comm <- n_comm ca_eig$pval <- pval# Store results simres$cor <-combine_cor(simres$cor, ressim) simres$eig <-rbind(simres$eig, ca_eig) }# Keep one set of simulated data simres$example[[parms$id]] <- sim } })# Save resultsaveRDS(simres, file = simu_files[[delta_name]])# Store result in list simres_list[[delta_name]] <- simresmessage(paste0("Output of system.time for delta = ", delta, ":"))print(rtime)message("The simulation took ", round(rtime[3]/60, 2)," minutes to run.") } else { # Use stored objectfor(delta in delta_list) { simres_list[[delta_name]] <-readRDS(simu_files[[delta_name]]) } }}
boxplots_list <-list()for(exp_plot inunique(dfsd_level$exp)) { dfsd_exp <- dfsd_level |>filter(exp == exp_plot) dfmean_exp <- dfmean_level |>filter(exp == exp_plot) legend_title <- legend_df$legend[legend_df$exp == exp_plot]# ---# Plot matrices# ---# Get ids to plot# We get the correspondence between ids and the parameter that changes fac_ids <-unique(dfsd_exp[, c("id", exp_plot)])# Sort to match the factor order sortfac <-order(factor(fac_ids[[exp_plot]])) ids <- fac_ids$id[sortfac]# Get corresponding data dats <- simres$example[ids]# Get color palette pal <-viridis(n =length(ids), begin =0, end =0.8, option ="D")names(pal) <- fac_ids[, 2]# Get plots widths/heights dims <-lapply(dats, function(x) dim(x$comm)) heights <-sapply(dims, function(x) x[1]) widths <-sapply(dims, function(x) x[2])# Get abundances and quantiles abds <-sapply(dats, function(x) max(x$comm[x$comm !=0])) cbounds <-sapply(dats, function(x) quantile(x$comm[x$comm !=0], probs =c(0.01, 0.99)))colnames(cbounds) <- ids gmatlist <-list()for (id in ids) {# Get corresponding data dat <- dats[[id]] comm <-as.data.frame(dat$comm)# Get the color colid <- pal[as.character(fac_ids[fac_ids$id == id, 2])]# Get the factor exp_fac <- fac_ids[fac_ids$id == id, exp_plot]# Perform multivariate analyses ca <-dudi.coa(comm, scannf =FALSE, nf =2) comm_reordered <- comm[order(ca$li[, 1]),order(ca$co[, 1])]# Censor data comm_reordered[comm_reordered >=max(cbounds[2, ])] <-max(cbounds[2, ]) comm_reordered[comm_reordered <=min(cbounds[1, ]) & comm_reordered !=0] <-min(cbounds[1, ]) gmat <-plot_matrix(comm_reordered, max_size = abds[id]/max(abds)*facsize,alpha = alpha_abund,col = colid) +theme(axis.text =element_blank(),axis.text.x.top =element_blank(),axis.title =element_blank(),axis.ticks =element_blank(),plot.margin =margin(t =0, r =10, b =0, l =10, unit ="pt"),panel.border =element_rect(colour = colid, linewidth =1),legend.position="none") gmatlist <-c(gmatlist, list(gmat)) }# ---# Plot boxplots# --- g1 <-ggplot(dfmean_exp) +ggtitle("Niche optima (CA ordination)") g2 <-ggplot(dfsd_exp) +ggtitle("Niche breadth (standard deviation)") facet_lab <-c("realized"="Realized niche", "fundamental"="Fundamental niche") boxplots <- g1 / g2 +plot_layout(guides ='collect') &ylim(c(0, 1)) &geom_boxplot(aes(y = cor, group =interaction(axis, type, .data[[exp_plot]]), x =factor(axis),col =factor(.data[[exp_plot]]))) &guides(colour =guide_legend(title.position ="top")) &scale_x_discrete(labels =c("1"="Axis 1", "2"="Axis 2")) &facet_grid(cols =vars(type), labeller =as_labeller(facet_lab)) &theme_linedraw() &ylab("Correlation (absolute value)") &theme(legend.position ='bottom',axis.title.x =element_blank(),legend.key.width =unit(0.05, "npc"),legend.key.spacing.x =unit(keyspacing[[exp_plot]], "npc"),legend.justification ="center",strip.background =element_rect(fill ="white"),strip.text =element_text(color ="black"))if (exp_plot =="abund_distri") {# Change labels boxplots <- boxplots &scale_color_manual(legend_title, labels = lab_abund_distri,values = pal) } else { boxplots <- boxplots &scale_color_manual(legend_title,values = pal) } matplots <-wrap_plots(gmatlist, nrow =1,heights = heights, widths = widths) boxplots <- boxplots /guide_area() /wrap_elements(matplots) +plot_layout(heights =c(3, 3, .5, 1.5)) boxplots_list[[exp_plot]] <- boxplots}
Code
for (n innames(boxplots_list)) {print(boxplots_list[[n]])ggsave(file.path(figures_path,paste0("re_fu_boxplots_delta", delta, "_", n, ".png")),width =15, height =18,units ="cm",dpi =600)}
References
Thioulouse, Jean, and Daniel Chessel. 1992. “A Method for Reciprocal Scaling of Species Tolerance and Sample Diversity.”Ecology 73 (2): 670–80. https://doi.org/10.2307/1940773.